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Abstract 

We compute the next-to-leading order QCD corrections to the production of two Z-bosons in 
the annihilation of two gluons at the LHC. Being enhanced by a large gluon flux, these corrections 
provide distinct and, potentially, the dominant part of the N^LO QCD contributions to Z-pair 
production in proton collisions. The gg —)• ZZ annihilation is a loop-induced process that receives 
the dominant contribution from loops of five light quarks, that are included in our computation in 
the massless approximation. We find that QCD corrections increase the gg ^ ZZ production cross 
section by 0(50% — 100%) depending on the values of the renormalization and factorization scales 
used in the leading order computation, and the collider energy. The large corrections to gg ^ ZZ 
channel increase the pp ^ ZZ cross section by about six to eight percent, exceeding the estimated 
theoretical uncertainty of the recent NNLO QCD calculation. 
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I. INTRODUCTION 


Production of pairs of vector bosons in proton collisions is one of the most interesting 
processes studied by ATLAS and CMS during the LHC Run I P-[3]. Indeed, pp —)■ ZZ, 
pp —)■ W^W~, and pp —)■ 77 were instrumental for the discovery of the Higgs boson. As 
the focus of Higgs physics shifts from the discovery to precision studies of the Higgs boson 
properties, di-boson production processes become essential for constraining anomalous Higgs 
boson couplings, for measuring the quantum numbers of the Higgs boson and for studying 
the Higgs boson width, see Refs. IlHZ]. Additionally, these processes provide important 
tests of our understanding of the Standard Model and can be used to constrain anomalous 
electroweak gauge boson couplings. 

Production of electroweak gauge boson pairs occurs mainly due to quark-antiquark annihila¬ 
tion qq —)■ ViV 2 - This contribution is known through next-to-next-to-leading order (NNLO) 
in perturbative QCD [SIRS]. However, as was pointed out in Refs. HMS], there is a sizable 
contribution from the gluon annihilation channel gg —)■ V 1 V 2 , whose signihcance depends 
on the selection cuts. For example, aggressive cuts applied to pp —)■ W~^W~ to separate 
the Higgs boson signal from the continuum background can increase the fraction of gluon 
fusion events in the background sample nzi- Since gg ^ V1V2 is a one-loop process and 
since production of electroweak boson pairs at leading order (LO) occurs only in the qq 
channel, the gluon fusion contribution to pp —)■ V 1 V 2 through NNLO only needs to be known 
at leading order, i.e. the one-loop approximation. Thus, all existing numerical estimates 
of the signihcance of the gluon fusion mechanism in weak boson pair production ignore ra¬ 
diative corrections to gg ^ ZZ that are, potentially, quite large ng. The need to have an 
accurate estimate of QCD corrections to gluon fusion processes for the Higgs width da ED] 
and generic off-shell measurements [2T1 - 123] was strongly emphasized in Ref. [7]. 

In this paper, we will focus on the calculation of the next-to-leading order (NLO) QCD 
corrections to the gluon fusion contribution to pp —)■ ZZ process. The largest contribution 
to gg ^ ZZ comes from quarks of the hrst two generations; these quarks can be taken to be 
massless. The situation is more complicated for quarks of the third generation. Ideally, we 
would like to include the (massless) bottom quark contribution and ignore the contribution 
of the massive top quark since, at leading order, the top-quark contributions change the 
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cross section by only about 1% (cf. Refs. [211 [25]).^ We can separate bottom and top 
contributions everywhere except in triangle diagrams that involve anomalous correlators of 
vector and axial currents. In these triangle diagrams, when bottom and top contributions are 
combined, the residual contributions are suppressed by the top quark mass, provided that we 
can assume it to be larger than any other energy scale in the problem. Unfortunately, in these 
diagrams top and bottom contributions can not be separated because the resulting theory 
is anomalous. To deal with this issue, we adopt the following strategy: we include quarks of 
the hrst two generations and the 6 -quark in our calculation in the massless approximation 
and we neglect all triangle diagrams whose contribution is then naturally associated with 
the quark contributions to gg ^ ZZ process. We note that the evaluation of the NLO QCD 
corrections to top quark mediated contribution to gg —)• ZZ process is not yet possible 
because the relevant two-loop amplitudes are not available. However, such contributions 
were recently studied in Ref. [26] in the approximation of a very large mass of the top quark. 
In that calculation quite large QCD corrections were found. 

Computing NLO QCD corrections to gg —)■ ZZ process is challenging because it is loop- 
induced. For this reason, the NLO QCD computation requires two-loop virtual matrix 
elements for gg ^ ZZ and one-loop matrix elements for gg —)■ ZZg processes. The recent 
progress in calculating two-loop integrals with two massless and two massive external lines 
[27H3I] made it possible to compute the required two-loop scattering amplitudes. Such 
amplitudes were calculated recently for qq —)■ V 1 V 2 [321133] and gg —)■ U 1 V 2 [331 [33] processes. 

The second ingredient that we need is the gg —>■ ZZg amplitude. Since this is a one- 
loop amplitude, it can be calculated in a relatively standard way, at least as a matter of 
principle. In fact, such calculations were performed in the past [3311371 and used to predict 
the production cross section ioi pp ^ ZZ + j. Automatic tools for one-loop computations 
can also deal with this process [381139] . Nevertheless, it is a non-trivial computation since, 
if we aim at calculating the NLO QCD corrections to gg ^ ZZ ^ 41, we require fast and 
stable calculation of helicity amplitudes for gg —ZZg process that includes decays of Z- 
bosons to leptons and can be extrapolated to soft and collinear kinematics of the final state 
gluon. Because of that, we decided to construct our own implementation of the scattering 

^ Contribution of the top quark loop becomes non-negligible in the region of high four-lepton invariant 
masses m^i > 2to(. 
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amplitude for gg —)■ ZZg using the unitarity methods 


The paper is organized as follows. In Section we present a brief review of the calculation 


of the two-loop scattering amplitude for gg ^ ZZ process. In Section III we discuss the cal¬ 
culation of the one-loop helicity amplitudes for gg —)■ ZZg and present numerical results for 


a kinematic point. In Section W we present numerical results for gg ^ ZZ contribution to 
pp ^ ZZ process at 8 and 13 TeV LHC at leading and next-to-leading order in perturbative 
QCD. We conclude in Section [V| 


II. THE TWO-LOOP SCATTERING AMPLITUDES FOR gg ZZ 

We start with a brief discussion of the two-loop scattering amplitudes for gg ^ ZZ process. 
Helicity amplitudes for this process were recently computed in Refs. [Ml [35] . In these 
references, each of the two independent helicity amplitudes for the process gg ^ ZZ ^ 41 
was written as linear combinations of nine form factors that depend on the Mandelstam 
invariants of the “prompt” process gg ^ ZZ and the invariant masses of the two Z bosons. 
The form factors are expressed in terms of polylogarithmic functions, including both ordinary 
and Goncharov polylogarithms. 

In this paper we use the results of Ref. |35] which are implemented in a C-I--I- code that 
can produce numerical results with arbitrary precision. In order to detect possible numer¬ 
ical instabilities, the code compares numerical evaluations obtained with different (double, 
quadruple and, if required, arbitrary) precision settings. If the results differ beyond a chosen 
tolerance, the precision is automatically increased. Of course, switching to arbitrary preci¬ 
sion increases the evaluation time substantially. Fortunately, we found that for phenomeno¬ 
logically relevant situations, the number of points where the code switches to arbitrary 
precision is negligible. Such points originate from kinematic regions where the two Z-bosons 
have either vanishing kinetic energies or vanishing transverse momenta. The amplitude 
squared is integrable in both of these regions, but, in practice, it can become numerically 
unstable. Since the contribution of these regions to the gg ^ ZZ cross section is relatively 
small, cutting them away, in principle, leads to an opportunity to perform stable numerical 

^ For recent reviews see Refs. |451 115] . 
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integration of the two-loop virtual correction over the four-lepton phase-space, resorting to 
quadruple precision only. However, we found that the improvement in performance achieved 
by cutting away the problematic regions is rather limited, so we used the default arbitrary 
precision implementation of the two-loop amplitude in practice. 

Since the gg ^ ZZ amplitude is one of the most complicated amplitudes that are currently 
known analytically, it is interesting to point out that the required evaluation times are 
acceptable for phenomenological needs. Indeed, calculation of all helicity amplitudes requires 
about two seconds per phase-space point in quadruple precision and, since the phase-space 
for gg ^ ZZ is relatively simple, one does not need excessively large number of points to 
sample it with good precision. 

For further reference we provide numerical results for the finite remainder of the one- and 
two-loop scattering amplitudes defined in gt-subtraction scheme, see Ref. |35]. The numerical 
results are presented for the choice of the renormalization scale /x = ^/s, where s is the 
partonic center-of-mass energy squared. The gj-subtraction scheme m is discussed in detail 
in Ref. |1H]. We consider the kinematical point 

g{pi) + g{P 2 ) (^/7)(P34) + (^/7)(P56) ^ e"(p 3 ) + e+(p4) + /i"(P5) + 


with (in GeV units) 


Pi = (99.5173068698129,99.5173068698129,0,0), 

P2 = (99.5173068698129,-99.5173068698129,0,0), 

P3 = (45.1400347869485,43.4878610174890, -9.85307698310431, 7.02463939683013), 

( 1 ) 

PA = (55.6586029753540, -27.4053916434553,48.1951275617684,4.90451560725290), 

P5 = (36.2015682945089, 34.5902512456859, -8.01242197258994, 7.06180995747356), 

Pa = (62.0344076828144, -50.6727206197196, -30.3296286060742, -18.9909649615566), 
and define a normalized amplitude through the following equation 


= WW X 10-» X ^ 67,67") 


Ai ,A2 ,Ae ,Aii 


dLIPS4. 


( 2 ) 


Note that in Eq.(|^ all the color factors have been factored out and dLIPS 4 is the standard 
Lorentz-invariant phase-space of the four final leptons. The color-stripped amplitude admits 
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Helicity amplitude 

1-loop 

2-loop 

A(l-,2-;3-,4+,5-,6+) 

-3.6020208 - 0.80680028 i 

-87.785548 -h 35.086257 i 

A(l-,2+;3-,4+,5-,6+) 

-h0.2507409 -h 0.38426042 i 

-hl8.585086 -h 7.5961902 i 


Table I: Results (in GeV for normalized qt remainder of one- and two-loop amplitudes for different 
choices of gluon and lepton helicities, evaluated at the scale /r = ^/s. See text for details. 


an expansion in the strong conpling constant 
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(3) 


Nnmerical resnlts for the two independent helicity amplitndes at one- and two-loops are 
given in Table [Tj We emphasize that the resnlts in Table |T] are given in the gt-snbtraction 
scheme, c.f. Ref. [35] . 


III. THE ONE-LOOP SCATTERING AMPLITUDE 0 ^ gggZZ 


In this Section, we discnss the compntation of the one-loop scattering amplitnde reqnired for 
the calcnlation of the inelastic process gg ^ ZZ + g.^ To this end, we consider the process 
0 9 {Pi) 9 {P 2 )g{P 3 )Z{p 45 )Z{pQj). Decays of the Z-bosons are allowed bnt, since we are 

interested in the on-shell prodnction of the two Z-bosons, we do not inclnde single resonant 
diagrams where one of the Z-bosons is emitted from the decay prodncts of the other one, 
see Fig. [Tj We will refer to the decay prodncts of the Z-boson with momentnm P 45 as the 
electron and the positron with momenta p 4 and ps and to the decay prodncts of the Z boson 
with momentnm p^j as the mnon and the anti-mnon with momenta pg and py, respectively. 
All leptons are taken to be massless. Since helicities of massless leptons are conserved, we 
only need to specify helicities of the hnal state leptons e~ and the allowed helicities of 
the positron and the anti-mnon in the hnal state are then antomatically hxed. 


^ To simplify the notation, in this section we do not consider photon-mediated four-lepton amplitudes. For 
phenomenological results discussed in Sec. IV we consider the full gg —>■ {Zj^){Zjg amplitude. 
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b) 

Figure 1: Representative Feynman diagrams for the 0 —)• gggZ{-^ e~e^)Z{-^ g,~g,'^) amplitude. 
Double resonant diagrams (a) are relevant for both the on-shell and the off-shell production. Single 
resonant diagrams (b) are only relevant for the off-shell production and are not included in our 
computation. See text for details. 


We write the interaction vertex of the Z-boson and a fermion pair as 

yf 7/.(1 + 75) , 7m(1-75) X 

^/7m/ e Qlj -^-+ Qrj -^- , / e (/, q). 

The left and right couplings for leptons and quarks are given by an identical formula 

_Vf±Af 

9l{R)j - ——^—, 
cos UT^ 


(4) 


(5) 


where we use i) Vi = —1/2-1- 2siv?6w, Ai = —1/2 for charged leptons; ii) 14 = 1/2 — 
4/3 sin^ 6w, A.'ii — 1/2 for up type quarks, and 14 — 1/2 2/3 sm Ad — 1/2 for 

down-type quarks. 

The 0 —)■ gggZZ scattering amplitude can be written as a sum of two terms 


= ghw (Tr 4f| + Tr Afi) , 


( 6 ) 


with Tr(f“ The two color-ordered amplitudes, stripped of their couplings to 

leptons and quarks, are defined as 


^ijk ~ (s'll Xj, Xk] Xe, A^) -|- gnuA^jf^ (A*, Xj, A^; Ae, A^)). 

In Eq.Q we introduced 

Cx^i = Dz{m]i) {gL,iSx- + 9r,iSx,+) , 


(7) 


( 8 ) 


where Dz{s) is the function related to the Breit-Wigner propagator Dz{s) = s/{s — -|- 

iMzTz). The couplings gff and g^^ are expressed through Z-boson couplings to quarks 
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propagating in the loops 


(9) 


dll = 

9 


ail = ^9i 

Q 


,1' 


Given these dehnitions, it is easy to see that the helicity amplitudes Xj, A^; Ae, A^) 

can be calculated for vector couplings of Z-bosons to leptons and quarks, provided that one 
keeps left-handed (right-handed) quarks propagating clockwise in the fermion loop. This is 
a natural separation if the scattering amplitudes are computed using the unitarity methods 
|40H44] ■ There is a useful relation between left- and right-handed helicity amplitudes for two 
orderings of external gluons 


-^ 3 , A2; Ae, A^) — —^f2^(Ai, A2, A 3 ; Ae, A^), 
-^ 3 ) A2; Ae, Xfj) = —Ai 23 {Xl, A2, A3; Ae, A^). 


( 10 ) 


These equations suggest that it is sufficient to compute LL and RR amplitudes for a single 
ordering; once this is done, all relevant amplitudes for the second ordering can be con¬ 
structed. Finally, we emphasize that we exclude the Breit-Wigner factor^ for the Z-bosons 
from the dehnition of the color-ordered helicity amplitudes but we include the 1 /s factor in 
its place; this can be clearly seen from the dehnition of the Dz{s) function in Eq.([^. 


It is well-known that any one-loop amplitude can be written as a linear combination of 
one-loop scalar integrals that include four-, three- and two-point functions and a rational 
part 


aLL^RR / \ \ \ . \ \ A _ \ LL,RR j . jy 

‘^ijk ‘^31 ^ G li it 


LL,RR 


( 11 ) 


The coefficients Cj in the above equation, as well as the rational part, can be calculated using 
unitarity methods. 


The idea of the unitarity method is that one can calculate the different discontinuities of the 
left- and right-hand sides of Eq. ( [IT| ) and then combine them in such a way that coefficients q 
are extracted algebraically. Calculation of the reduction coefficients and the rational part can 
be performed either analytically or numerically. In this paper, we use a mixed approach. 
We compute the coefficients c* using numerical four-dimensional unitarity introduced in 


As we mentioned earlier, we are interested in the on-shell production of the two Z-bosons in this paper. 
However, we construct the relevant piece of gg —>■ ZZ amplitude in full generality, including Breit-Wigner 
propagators for the Z-bosons, to enable its later use to study QCD corrections to gg —t ZZ* process. 



Helicity 

LL 

RR 

-4i23(+) +) +; —) —) 

-42.714233 -h 117.60020 i 

-138.32358-h 139.68765 i 

-4i23(+) +) —; —) —) 

-M34.26016-h 161.13392 i 

-^138.09750-h 188.27580 i 

-4i23( —) —; +; —) +) 

-32.287418-^2.1139258 i 

-31.55258-F 32.433444 i 


Table II: Results (in GeV for normalized color-ordered amplitudes for 0 —?■ gggZ{e~^e g ) 

process, for different choices of gluon and lepton helicities. See text for details. 

Ref. |13]. The rational part, on the other hand, is computed analytically using the method 
described in Refs. 150] . Technical details about the unitarity methods used for one-loop 
computations in QCD can be found in Ref. 0SI- 

From the point of view of the unitarity methods, the peculiarity of 0 —)■ gggZZ process is 
that it involves two colorless particles, making full color ordering for scattering amplitudes 
impossible. This has the following implications. Any unitarity computation starts with 
the list of independent “parent diagrams” that are subsequently cut into on-shell scattering 
amplitudes. Although parent diagrams are independent by construction, not all their cuts 
are, if permutations of external particles are allowed. The challenge, therefore, is to start 
with the “parent diagrams”, write down all the cuts that they might have and then exclude 
all the cuts that are not independent. This issue was successfully dealt with in the context 
of many recent calculations of one-loop scattering amplitudes for quarks, gluons and vector 
bosons, see e.g. Refs. In this paper, we construct the independent set of unitarity 

cuts following the methodology explained in those references. 

After identifying independent cuts, we hnd 39 quadruple, 45 triple and 18 double cuts. 
There are no single-line cuts since internal fermions in our calculation are massless. Each 
of these cuts is described by a product of tree-level color-ordered amplitudes. The required 
helicity amplitudes include qgq, qZq, qggq, qgZq, qZZq, qgggq, qggZq, qgZZq and qgggZq. 
Here, we use a generic notion of a Z-boson for an external vector particle but what we really 
mean are amplitudes with the vector current sandwiched between lepton and anti-lepton 
spinors. The relevant tree-level amplitudes can be extracted from different publications; we 
have mostly benehted from a comprehensive description of helicity amplitudes that involve 
quarks, gluons and vector bosons in Ref. |55) . 
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The calculation of the rational part of the 0 —?• gggZZ amplitude is performed analytically, 
using techniques suggested in Ref. jlHlEn]- Similar to the cut-constructable part the rational 
amplitude receives contributions from quadruple, triple and double cuts. However, for the 
case of gg —)■ ZZg amplitudes, the double cut contribution vanishes; the rational part 
therefore can be reconstructed from the calculation of boxes and triangles. Unfortunately, 
even in this case, the analytic results for the rational part are unwieldy and we choose not 
to present them here. 


For further reference, we give numerical results for the scattering amplitudes below. We 
consider a kinematic point (momenta are given in GeV) 


Pi = (-238.714576090637,-238.714576090637,0,0), 

P2 = (-1021.22119318758,1021.22119318758,0,0), 

P3 = (250.736037681104, -207.896850811885, -124.613643938661,64.1786550096635), 
Pa = (553.889863453468, -495.644737899924, -245.099246845329,32.5059044554765), 
P5 = (91.0664644166627,49.0057636944973, 76.1125676676337, -9.92033815503652), 
P6 = (197.326337775966, -3.11006048502754,183.877222508616, -71.5344542606618), 

P7 = (166.917065951017, -124.860731594604,109.723100607740, -15.2297670494417), 
and define a normalized primitive amplitude through the following equation 


( 12 ) 


A 


LL,RR( 


ijk (^i. K, Ap) = X 10 " X xt,?;.t®"(A,, Aj, \f, K, A„). 


(13) 


The results for certain helicity combinations of gluons and leptons are given in Table |TT| 
We emphasize that diagrams where one Z-boson is emitted by decay products of another 
Z-boson, see Fig. [^b), are not included in our calculation. The result for the amplitude 
squared and summed over colors and helicities of gluons and leptons was checked against 
the results of the OpenLoops program [2B] for a large number of kinematic points.^ Finally, 
we note that the evaluation of the amplitude squared, summed over color and helicities, 
takes about 0.1 seconds per phase-space point, making our implementation adequate for 
phenomenological needs. 


In the context of NLO QCD computations, the process gg —>■ ZZ + g represents an inelastic 
contribution. This inelastic contribution should be integrated over all energies and angles of 


^ We are indebted to J. Lindert for making this comparison possible. 
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rricut [GeV] 7714/ [GeV] 

Figure 2: Up, left: cumulative cross section for gg —)• [Zjj —>■ e^e~g'^g~ at the 8 TeV 
LHC as a function of the lower cut on four-lepton invariant mass. Up, right: distribution of the 
invariant mass of the four leptons in the reaction gg —)• [Zj^){Zj —)> at the 8 TeV 

LHC. Lower panes show ratios of the LO (yellow) and NLO (blue) distributions evaluated at three 
different scales to the LO distribution evaluated at g = 2mz- Low: same as above for the 13 TeV 
LHC. 

the emitted gluons, including the vanishingly small ones. Calculation of one-loop amplitudes 
for gg —)■ ZZg process becomes unstable if the gluon in the final state becomes soft or 
collinear to the collision axis. We deal with these instabilities by switching to quadruple 
precision where appropriate. To obtain the gg —)■ ZZ cross section through NLO QCD, 
we combine elastic and inelastic contributions using the g^subtraction m and, as a cross¬ 
check, the FKS subtraction [56] methods. The results that we present in the next Section 
are obtained by combining computations performed using the two subtraction schemes. 
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IV. NUMERICAL RESULTS 


In this Section we present the results of the calculation. We consider the process gg —)■ 
{Z/'y){Z/'-j) —>■ at the LHC.® We generate invariant masses of Z bosons around 

rriz, using Breit-Wigner distributions. We require the e^e“ and gZgT pair to have invariant 
masses mu G (60,120) GeV. We use leading (next-to-leading) order parton distribution 
functions and the strong coupling constant for one- and two-loop calculations, respectively. 
We employ the NNPDF3.0 set of parton distribution functions and obtain the relevant values 
of the strong coupling constant from NNPDF routines |58) . 

We begin with presenting the results for the total cross sections at the 8 TeV LHC. We hnd 

= 0.97«| fb, = 1.8«| fb, (H) 


where the central values refer to the renormalization and factorization scales set to g = 2mz 


and the upper (lower) values to g = mz {g = 4m^). It follows from Eq.(14) that QCD 
corrections to gg ^ ZZ are large - the NLO cross section increases the LO cross section by 
0(60% — 110%), depending on the renormalization scale. For g = 2mz, the cross section 
increases by 85%. 


A similar situation occurs at the 13 TeV LHC. We hnd 


= 2.81“;' fb, = 4.7«:1 fb. (15) 

The NLO QCD corrections to gg ^ ZZ at 13 TeV LHC are again signihcant but somewhat 
smaller than corrections at 8 TeV. Indeed, for the central value of the renormalization and 
factorization scales g = 2m^, the cross section increases by 67%. For other values of the 
renormalization and factorization scales, the cross section increases by 0(40% — 90%). 

The large size of the QCD corrections is reminiscent of the large QCD corrections to Higgs 
production in gluon fusion gg ^ H [59]. In addition, similar to the Higgs production 
case, the scale variation of the leading order cross section provides a poor estimate of the 
magnitude of next-to-leading order corrections |59|. We note that if we take proximity of 

® We remind the reader that we only include double resonant diagrams, see Fig|^a). Single resonant 
diagrams Figj^b) are only relevant for far off-shell production. They can be obtained by appropriate 
modifications of the gg Zg amplitudes, see e.g. m- 
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radiative effects in gg — )■ ZZ and gg ^ H seriously, we should probably take /i = (2m^)/2 = 
rriz as the scale for which higher-order radiative corrections to gg ^ ZZ will most likely be 
small. Thus, our best estimates ot gg —)■ ZZ contributions to pp ^ ZZ production cross 
section at 8 TeV and 13 TeV LHC are 

'^pp^ZZ (8 TeV) = 2.0(2)fb, TeV) = 5.1(4)fb. (16) 


Our results have important implications for the recently computed NNLO QCD corrections 
to pp —)■ ZZ [TUI [TU] at the 8 TeV LHC. In that case, the NNLO QCD corrections computed 
at the scale p = mz turned out to be close to 15%. However, a signihcant fraction - 60% 
of the total NNLO QCD correction - is due to the leading order contribution gg — ZZ. 
Our current computation shows that gg ^ ZZ receives large radiative corrections and the 
natural question is how these hndings affect the central value oi pp —)■ ZZ cross section 
obtained in Refs. dQiEg and the theory uncertainty assigned to it. 


To answer this question, we note that in Refs. [ini[i2] the central scale was chosen to be 
p = mz and that NNLO parton distribution functions were used for the calculation of 
gg ^ ZZ cross section. Relative to our choices, the lower renormalization and factorization 
scale increases the cross section while the choice of NNLO parton distribution functions 
makes the cross section smaller. We re-computed the LO gg ^ ZZ cross section using the 


setup of Ref. im and compared it with our best value given in Eq.([l6|. We hnd that, to 
match our best prediction, the 8 TeV gg ^ ZZ cross section of Ref. ra should be increased 
by about 80%. In turn, this will lead to an increase in the total NNLO QCD correction to 
pp ^ ZZ aX ^ TeV from the current 12%, as calculated in Ref. [TU], to 18%. This increase is 
beyond the 0(3%) scale variation of the NNLO QCD result for pp —>■ ZZ used in Ref. [I0| 
to estimate the current uncertainty in the theoretical prediction for pp ^ ZZ cross section. 
Similar arguments also apply at the 13 TeV LHC. In this case the 16% corrections quoted 
in Ref. [1^ would increase to approximately 23%. 


Next, we consider kinematic distributions. We begin with the invariant mass distribution 
of the four leptons produced in gg ^ ZZ shown in Fig. While radiative corrections are 
significant for all values of 777.4;, they become smaller at higher values of four-lepton invariant 
masses. This is clearly seen in Fig. for both differential and cumulative^ cross sections 


^ For different cuts on 7714;. 
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Figure 3: Left: transverse momentum distribution of an e^e“ pair at the 13 TeV LHC. Right: the 
hardest lepton transverse momentum distribution at the 8 TeV LHC. Lower panes show ratios of the 
LO (yellow) and NLO (blue) distributions evaluated at three different scales to the LO distribution 
evaluated at /U = 2mz- 

and for both the 8 TeV and the 13 TeV LHC. This result is important for studies of the 
Higgs off-shell production where good understanding of the shape of four-lepton invariant 
mass distribution is an important pre-requisite for constraining the Higgs width. Note that 
for rriii > 2mt top-quark contributions, neglected in our computation, become relevant. 

In Fig.j^we show the transverse momentum distributions of the pair and of the hardest 
lepton in the event. The QCD corrections to the transverse momentum distribution of the 
e’''e“ pair decrease for large values of p±^e+e-i similar to what is seen in the four-lepton 
invariant mass distribution. On the other hand, the QCD corrections for the transverse 
momentum distribution of the hardest lepton are independent of the lepton p±. 


V. CONCLUSIONS 

In this paper we computed QCD corrections to the production of a pair of Z-bosons in 
gluon fusion through loops of massless quarks. We found that QCD corrections are large; 
they change the production cross section by almost a factor of two. These large QCD 
corrections are in line with expectations that transition of two gluons to a colorless hnal 
state is strongly affected by QCD radiative effects; QCD corrections of similar magnitude 
were observed earlier in theoretical calculations of gg ^ H [3^] and gg —?• 77 110] cross 
sections. 
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Large QCD corrections to gg ^ ZZ are important for a nnmber of reasons. First, since 
the gg ^ ZZ process provides a significant fraction of the NNLO QCD contribntion to 
pp —)■ ZZ, onr resnlt snggests that existing theoretical predictions for pp ^ ZZ shonld 
be increased by six to eight percent, depending on collider energy. Since snch an increase 
in the central valne is ontside the existing estimates of the residnal theory nncertainty of 
pp ^ ZZ cross section, it will have important conseqnences for ongoing comparisons of 
experimental and theoretical resnlts for pp ^ ZZ at the LHC. Second, good nnderstanding 
oi gg ^ ZZ at high fonr-lepton invariant masses is crncial for the so-called off-shell stndies 
of the Higgs boson and, in particnlar, for the indirect determination of its width. The 
NLO QCD calcnlation of gg ^ ZZ process allows ns to predict the gg ^ ZZ contribntion 
to pp —)■ ZZ cross section and kinematic distribntion with the precision of abont 10%; 
this implies a residnal theoretical nncertainty on pp —)■ ZZ cross section of jnst abont two 
percent. Snch a small nncertainty in the fonr-lepton prodnction cross section is an essential 
prereqnisite for the snccess of forthcoming off-shell stndies of the Higgs boson, see a related 
discnssion in Ref. |7j. 

As a hnal comment, we note that onr calcnlation opens np a nnmber of fntnre research 
directions. Indeed, it is interesting to extend onr calcnlation by combining massless and 
massive loop contribntions to gg ^ ZZ and by inclnding single resonant contribntions and 
the interference of prompt gg ^ ZZ and gg —?■ H* ZZ amplitndes. This will allow 
ns to explore the region of fonr-lepton invariant masses both below the threshold of ZZ 
prodnction and at very high invariant masses. We plan to do this in the near fntnre. 
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